Semiclassical Correlation in Time-Dependent Density Matrix Functional Theory 
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Lack of memory (locality in time) is a major limitation of almost all present time-dependent 
density functional approximations. By using semiclassical dynamics to compute correlation effects 
within a density-matrix functional approach, we incorporate memory, including initial-state de- 
pendence, as well as changing occupation numbers, and predict more observables in strong-field 
applications. 
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The impact of time-dependent density functional the- 
ory (TDDFT) [l|, [2| on calculations of excitation spec- 
tra and response in atoms, molecules, and solids is evi- 
dent in its increasing use. In such applications a weak 
perturbation is applied to the system beginning in its 
ground state, and usually the exchange-correlation (xc) 
effects are treated with a ground-state approximation. 
Generally the results are usefully accurate, but specific 
cases (e.g. charge transfer excitations, optical response of 
solids, etc) require improved approximations undergoing 
intense research. 

In principle, TDDFT also applies beyond the linear 
response regime, but success has been slower. There 
are three main reasons. First, many of the observables 
of interest are not simply related to the time-dependent 
one-body density: in addition to the approximation for 
the xc functional, new approximate "observable func- 
tionals" arc needed to extract the quantity of interest 
from the Kohn-Sham (KS) system. Even with an exact 
xc potential, they would remain elusive. These include 
double-ionization probabilities and momentum-densities, 
and naive approximations to these generally fail (y, |7J. 
Second, lack of memory dependence in the usual xc ap- 
proximations has been shown to be often far more prob- 
lematic than in the linear-response regime [3|, Q[. The 
exact functional depends on the history of the density as 
well as on the initial state. Different initial states lead to 
fundamentally different xc potentials i4] . But no approx- 
imation today has initial-state dependence, almost all ne- 
glect history-dependence, and all violate an exact condi- 
tion on memory-dependence, derived in Ref. [5]. Third, 
a particularly severe difficulty is encountered when a 
system starting in a wavefunction dominated by a sin- 
gle Slater determinant (SSD), evolves to a state that 
fundamentally needs at least two SSDs to describe it. 
This is the time-dependent (TD) analog of ground-state 
static correlation, and arises in electronic quantum con- 
trol problems [f| Q, in ionization [7), and in coupled 
electron- ion dynamics Q. The TD KS system evolves 
the occupied orbitals under a one-body Hamiltonian, re- 
maining in an SSD: the KS one-body density matrix is 
always idempotent (even with exact functionals), while, 
in contrast, that of the true system develops eigenvalues 
(natural occupation numbers) far from or 1 in these ap- 



plications. The exact xc potential and observable func- 
tionals consequently develop complicated structure that 
is difficult to capture in approximations. For example, in 
Ref. 7], a simple model of ionization in two-electron sys- 
tems showed that the momentum distribution computed 
directly from the exact KS system contains spurious os- 
cillations due to using a single, necessarily delocalized or- 
bital, a non-classical description of the essentially classi- 
cal two-electron dynamics. Ref. [5J discussed the unusual 
and non-intuitive xc potentials that arise in certain elec- 
tronic quantum-control problems, e.g. He Is 2 — > ls2p. If 
the overlap between the initial and final states is target- 
ted, the maximum that can be achieved is 0.5 |9(], while 
close to 0.98 is achieved for the true interacting problem. 

Recent pioneering strides in TD density-matrix func- 
tional theory (TDDMFT) show this alternative approach 
can overcome some of the challenges of adiabatic TDDFT 
in linear response [10j . e.g. adiabatic TDDMFT func- 
tionals were shown to capture charge-transfer excita- 
tions well. All one-body observables are directly ob- 
tained. However, adiabatic functionals bootstrapped 
from the usual ground-state DMFT disappointingly can- 
not change occupation numbers unless some unusual 
structural changes are made in the form of the func- 
tional [10]. The first real-time TDDMFT calculations 
have been performed recently [ill ] , using an extra energy- 
minimizing procedure at each time that results in time- 
dependent occupation numbers. 

In this Letter, we present a new approach to correlation 
in TDDMFT that makes a significant step in solving all 
the problems mentioned above. We work within real-time 
TDDMFT and use a semiclassical approximation for the 
correlation term in the equation of motion while evaluat- 
ing the other terms exactly quantum-mechanically. All 
one-body observables are obtained directly, with correla- 
tion effects treated semiclassically. It is the first density- 
matrix(or density-) functional approach that has initial- 
state dependence, with memory naturally carried along 
by the classical trajectories, and the first real-time ap- 
proach that can change occupation numbers significantly 
away from the adiabatic limit. A heirarchy of semiclassi- 
cal approximations for the correlation term is discussed, 
decreasing in accuracy but also in computational cost. 
On the first level, correlation is obtained exactly to 0(h), 



while at the lowest level quantum mechanics enters only 
in the determination of the initial state, with the dy- 
namical correlation obtained via pure classical evolution. 
Despite its simplicity, we demonstrate via a simple ex- 
ample that this latter approach yields time-dependent 
occupation numbers. 

Semiclassics lies at the very foundations of the 
earliest density functional theories that predate the 
rigorous DFT of Refs. [12]. Its semiclassical origins 
were however somewhat forgotten in the developments 
in the 1990's of ground-state functionals, and only very 
recently have been reawakened [13| . Until now, semi- 
classical methods have not been applied to functional 
development in TDDFT nor in TDDMFT, although 
mean-field semiclassical methods have been used to 
approximate KS dynamics (e.g. Vlasov approaches to 
metal clusters in strong fields [l4J|). TDDMFT may 
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equally be viewed as a "phase-space-density functional 
theory", as reflected for example in the relation be- 
tween the one-body Wigner function w(r, p,i) and 
the spin-summed one-body density matrix pi(r',r, t) — 
N J2 n ..a N I d 3 r 2 ..d 3 r N ^*(r'a u x 2 ..x N ,t)^(ra u x 2 ..x N} t): 



w(r,p,t) = (J^j y'd 3 yp 1 (r- y /2,r+y/2,i) e i Py (1) 



(Here Xi = (r^, er^) indicates spatial and spin coordinates, 
and atomic units are used throughout). This observa- 
tion suggests the utility of semiclassical approaches, as 
we shall see shortly. 

All the one-body terms in the equation of motion for p\ 
can be treated exactly, and for spin-unpolarized systems: 



Vi(r',r,t) = (-V 2 /2 + v(r,t)+V 2 /2-v(r',t))p 1 (r\r 7 t) 

+ / d 3 r 2 / e(r,r',r 2 ) f n(r 2 ,t)pi(r',r,t) - ~pi(r' ) r 2) t)pi(r 2 ,r,t) J + / d 3 r 2 / oc (r,^r 2 )p 2c (r',r 2 ,r,r 2 ( 
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where we have decomposed the second-order density ma- 
trix, p 2 (r',r 2 ,r,r 2) i) = N(N ~ 1) £ CTl .. ffJV J d 3 r 3 ..d 3 r N 
**(r'cri, r 2 CT 2 , x 3 ..x N , t)^(rcri, r 2 a 2 , x 3 ..x N , t) 
pi(r[,ri,t)pi(r' 2 ,r 2 ,t) - p 1 (r' 1 ,r 2 ,t)p 1 (r 2 ,r 1 ,t)/2 + 
/ o 2c (r' 1 ,r 2 ,ri,r 2 ,t): the first term is the non-interacting, 
uncorrelated product, the second term takes care 
of the Pauli principle at the uncorrelated level, and 
the third term is the correlation component, whose 
functional dependence on p\ is unknown. In TD- 
DMFT, p 2 (r',r 2 ,r,r 2 ,i) is to be approximated as a 
functional of p\ and the initial interacting state "3/oj 
92 \p\ , ^o] ■ We have defined the electron-interaction 
kernel / oc (r, r', r 2 ) = 1/ |r — r 2 | — l/|r' — r 2 |. If we evolve 
the A-electron interacting system in a TD external 
potential v(r,t), that there is a one-to-one mapping 
between v(r,t) and p\ (or w(r,p,t)) for a given initial- 
state x I'o( r i--- r Jv) follows directly from the Runge-Gross 
theorem. The mapping holds also for external vector 
potentials [7| but we will focus on scalar potentials 
at present. An immediate advantage of replacing the 
coordinate-space density with the density-matrix as 
basic variable is that it directly gives the expectation 
value of any one-body operator: no additional observable 
functionals are needed for momentum-distributions or 
kinetic energies, for example. There is no KS equivalent: 
because of idempotency of non-interacting density 
matrices, it is impossible for a non-interacting system 
of electrons to have the same phase-space density as a 
system of interacting electrons. 

Ideally, the approximation made for p 2 [pi, ^o] captures 
correlation, memory-dependence including initial-state 



dependence, and, most importantly for the quantum 
control and ionization applications mentioned earlier, 
can yield time-dependent occupation numbers, /i(i), de- 
fined via the natural orbital decomposition pi(r, r',i) = 

The TDDMFT developments have so far been pre- 
dominantly within linear response [l(| , investigating adi- 
abatic functionals for p 2 . Our approach computes the 
correlation component of p 2 via semiclassical dynamics, 
focussing on full dynamics (not linear response). We 
propagate Eq. ^ treating all terms except the last ex- 
actly quantum-mechanically. The last term is treated as 
a driving term: we approximate p 2c s=s p 2 °, evaluated 
separately via semiclassical dynamics, calculated from 
running classical trajectories in the A-body interacting 
phase-space. 

Semiclassical methods construct an approximate quan- 
tum propagator utilizing classical trajectory information 
alone. Although there are a variety of forms, the essential 
structure is a sum over classical trajectories: 



cl.traj 



Ci(t)e 



iS,(t)/h 



(3) 



where Si(t) is the classical action along the ith trajec- 
tory, and the prefactor C'i (t) captures fluctuations around 
the classical path. Semiclassical approaches are able 
to capture quantum effects such as interference, zero- 
point energy, tunneling (to some extent), while gener- 
ally scaling favorably with the number of degrees of free- 
dom. Based on classical trajectories, intuition about 



the physical mechanisms underlying the dynamics can be 
gained. Although mostly applied to nuclear dynamics in 
molecules, there have been applications to electrons [15j . 
Scmiclassical formulae have been derived both from 
largely intuitive arguments (e.g. Ref.|16|h as well as 
from careful rigorous asymptotic analyses of the quan- 
tum propagator (see e.g. Refs.jl71|') that satisfy TDSE 
to order h. Miller [18[ showed the equivalence of differ- 
ent semiclassical representations within stationary-phase 
evaluation of the transformations. The most popular 
is the Heller-Herman-Kluk-Kay [H [H H3|, which is 
a "semiclassical rigorization" of Heller's frozen Gaus- 
sian approach, uniformly solving the TDSE to 0(h). It 
is a sum over initial points in (TV-body) phase-space, 



io 



^o'E ) 



((ri(0) j p 1 (0))...(r*(0),ptf(0))): 



(e-***) 



sc 



j2A/ r 



=0 



(27T) 



M 



|« t >Ct<g )e«-^< So 



where: z = (r(t),p(t)) obeys Hamilton's equations 
r = p(t), p = -V#(r,p,i), 



(4) 



(5) 



M = 3iV is the dimensionality of configuration space, 
and S(z , t) is the classical action, J (T — V)dt, for 
a trajectory which begins at the phase space point 

(r , p ), reaching point (r , p ) at time t. The state 
—0 =o — * =t 

|z) is a product of coherent states for each coordi- 
nate, labelled by their centers in phase-space: (x\q,p) = 

(^) exp(— j(x — q) 2 /(2) + ip(x — q)), where 7 is a 
chosen width parameter. The pre-exponential determi- 
nant C t (z) accounts for fluctuations about the classi- 
cal paths: when 7 is chosen identical for all particles, 
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Typically, the phase-space integral is performed via 
Monte-Carlo, with initial conditions weighted by the ini- 
tial wavepacket (z |^o) ■ Due to the evaluation of the 
prefactor C, the numerical effort per trajectory scales as 
N 3 ; methods which neglect this scale more favorably as 
N but at the cost of losing accuracy and semiclassical 
rigor. While Monte-Carlo integration scales as y/N for 
positive integrands, the phase-space integral can be dif- 
ficult to converge due to its oscillatory nature, especially 
for many degrees of freedom and chaotic dynamics, and 
so various sophisticated integral-filtering techniques, or 
"forward-backward" (FB) methods [21| have been for- 
mulated, allowing rigorous semiclassical calculations for 
up to 100 degrees of freedom |22j |. 

Applying Eq. (J4J to propagate ^/, then computing p 2 c 
via integration constitutes our highest level of semiclas- 
sical heirarchy for correlation. We compute 



where n sc (r 2 , t) is the one-body density obtained from 
the semiclassical calculation, the diagonal of Eq. (|6]). Fi- 
nally, this pic is input into Eq. J2} as a driving term. 
At this level, correlation effects are exact to 0(fi), while 
all other effects are quantum-mechanically exact. Diffi- 
culties with convergence of the highly-oscillatory integral 
in Eq. (J4]), and the iV 3 scaling of the pre- factor, render 
this impractical for more than a few electrons, to the 
point where, for many cases, little computational benefit 
is gained over running the full quantum mechanics. 

Instead the FB idea lends itself particularly well to 
our purposes: we can take advantage of significant phase- 
cancellation between the propagation of ^* and that of '3/ 
in calculating p 2 c . Applying the semiclassical propagator 
Eq. ([JJ to both the \& and ^>* appearing in p%, and doing 
some intermediate integrations via stationary phase, the 
second level in our semiclassical heirarchy is: 



pf (r',r 2 ,r,r 2 ) 
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(2tt) 3 



-j3 J j3 „/ 



3JV+2 / « lo" Z U a 



J(S(t)-S'(t)) 



x g(r',r, r 2 ; Z i, t z 2it z lit z 2 , t )(*oli )(i l*o) (8) 



Io 



and 



where 

(qi,t(-*)i Pi.t(-t): q 2 ,t(- f ). P2,< (-*)> r 3,o, p 3 ,o---rjv,o, PiV.o) 
0(r',r,r 2 ;zi it z^ t zi, t z 2l t) = 

(r / |zi jt )*(r2K )t )*(r|z 1)t }(r 2 |z 2 , t }. That is, an ini- 
tial phase-space point q , p is classically evolved for 

=0 =0 
time t, when the phase-space points of the the first 

two particles are shifted to (qj 1; pj 1; c^ t2 ,p' t2 ), before all 
particles evolve back to time zero. There is therefore 
significant cancellation of phase (S(t) — S'(t)), that 
would generally result in good convergence of Monte 
Carlo evaluation of this phase-space integral, even for 
many electrons. Further, the product of the numerically 
expensive pre-factors has been reasonably approximated 
to 1 for many electrons. The true initial state appearing 
in Eq. (jSJ) is in practise approximated by a few (KS) 
SSD's, or by a high-level wavefunction calculation, if 
a stationary state. Eqs. ^ and (J7J are then used to 
extract the semiclassical correlation component p 2 £, 
capturing interference and zero-point energy effects, that 
is then input into Eq. ([2]) as a driving term. 

An even more simple prescription is obtained by ne- 
glecting the phase and prefactor a ltog ether: this yields 
the quasiclassical Wigner method [23|], and can also be 
shown to result from a linearization of the FB 1211 : 
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! (|,p,t) = w N [L_ t ,g_ t ,t = 



(9) 



pf(r',r,i) = j^jJpT (r',r 2 ,r,r 2 ,t)d 3 r 2 



and extract 



p^(r',r 2 ,r,r 2 ,i) = pf (r', r 2 , r, r 2 

+ pf(r',r 2 ,i)pf(r 2 ,r,i)/2 



from which, by integration, a quasiclassical approxima- 
tion to the correlation component of p 2 is obtained, and 
inserted as a driving term into Eq. [5] This represents 
(q\ the lowest level of our semiclassical heirarchy: in com- 
puting the correlation, while scaling classically with the 
system size, all interference is lost, quantum mechanics 
enters only in determining the initial Wigner function, 
*) ~ Pi ( r 7 r > t) n ( r 2: ''and when the wavefunction becomes delocalized, this ap- 
proximation degrades. Nevertheless quasiclassical meth- 



ods (even of the entire dynamics) have proven useful in 
analyzing electron ionization distributions [241 ]. 

Our prescription thus results in a semiclassical approx- 
imation for the correlation component to pi in the equa- 
tion of motion Eq. (|2J for p\ , all other terms of which are 
treated exactly quantum-mechanically. But can our ap- 
proach lead to time-dependent occupation numbers? To 
illustrate this, we consider a simple model system, the 
two-electron Moshinsky atom [251 ] : 

H = - l -(Vl + Vt)+'^(rl + rl) + \(r 1 -r 2 f (10) 

Although a poor model of a real atom, its purpose here 
is simply to demonstrate that even the lowest level qua- 
siclassical approximation to correlation is able to cap- 
ture changing occupation numbers. Its harmonic nature 
makes it exactly solvable, and we apply a simple sinu- 
soidal force constant, k(t) = 1 — 0.05sin(2i), that encour- 
ages population transfer to the first accessible excited 
state (an excitation in the center of mass coordinate), 
from the initial ground-state, a spin-singlet. Moreover, 
due to its harmonic nature, the quasiclassical and semi- 
classical propagations are exact. Figure Q] plots the oc- 
cupation numbers, fi(t), of the spatial natural orbitals, 
obtained from diagonalizing pi(r, r', t) at each time t. In 
striking contrast, the results of the usual adiabatic ap- 
proximations in TDDMFT would yield constant occupa- 
tion numbers. How well our approach works for more 
realistic systems is currently being tested; this exam- 
ple however illustrates that it certainly does not suffer 
from the inability to change occupation numbers. Aside 
from the significance in quantum control problems, lack 
of time-dependent occupation numbers impacts observ- 
ables, e.g. the momentum distributions are qualitatively 
incorrect pj- 



FIG. 1: Occupation numbers for the model system: quasi- 
classical correlation shown is exact, while the usual adiabatic 
TDDMFT approximations yield constant straight lines. 




In summary, we have presented a semiclassical ap- 
proach to correlation in TDDMFT, that (i) naturally 
captures history-dependence and initial-state depen- 
dence (for the first time) at the semiclassical level, as 
memory is carried along with the classical trajectories 
composing p^ , (ii) directly yields all one-body observ- 
ables, and (iii) changes occupation numbers. Correlation 
is treated semiclassically, while all other terms determin- 
ing the density- matrix are exact. The highest semiclas- 
sical level yields correlation exactly to 0(h), but will be 
impractical in many cases of interest; the approximate 
semiclassical treatment (Eq. [8]) will still capture quan- 
tum many-body effects in a computationally efficient 
way. The simplest approximation (Eq. [9]) scales classi- 
cally, so is well worth investigating, especially since the 
other terms in the equation of motion for p\ are treated 
exactly. As there is no guarantee of iV-representability of 
p\ being preserved in the semiclassical dynamics (at least 
beyond 0(h)), tests on more realistic systems than that 
presented here are necessary. Treating many of the chal- 
lenging aspects of approximate density-functional meth- 
ods described earlier, it is a promising approach to study 
electron dynamics in strong fields. 
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